subroutine set_source_pos(gid,rank,lid,px,py,pz)
  use param
  use global_param,only:isqrt3
  use geom_data_param,only:geom_obj
  use label_data_param,only:obj_type
  implicit none

  integer,intent(in) :: gid,rank,lid
  real(dbl),intent(out) :: px,py,pz
  integer :: dir
  real(dbl) :: rn,rn2,d1,dist,dr

  dist=geom_obj(gid)%diss_rd(rank)
  dr=dist*isqrt3

  select case(geom_obj(gid)%geom_type)
  case(1)
     ! point
     px=geom_obj(gid)%x1
     py=geom_obj(gid)%y1
     pz=geom_obj(gid)%z1

     if (obj_type(lid)%does_move) then

        select case(obj_type(lid)%motion_type)
        case(0)
        case(1)
           call random_number(rn)
           dir=ceiling(rn*8.0)
           select case(dir)
           case(1)
              px=px+dr
              py=py+dr
              pz=pz+dr
           case(2)
              px=px-dr
              py=py+dr
              pz=pz+dr
           case(3)
              px=px-dr
              py=py-dr
              pz=pz+dr
           case(4)
              px=px+dr
              py=py-dr
              pz=pz+dr
           case(5)
              px=px+dr
              py=py+dr
              pz=pz-dr
           case(6)
              px=px-dr
              py=py+dr
              pz=pz-dr
           case(7)
              px=px-dr
              py=py-dr
              pz=pz-dr
           case(8)
              px=px+dr
              py=py-dr
              pz=pz-dr
           end select           
        case(2)
        case(3)
           call random_number(rn)
           dir=ceiling(rn*8.0)
           select case(dir)
           case(1)
              px=px+dr
              py=py+dr
              pz=pz+dr
           case(2)
              px=px-dr
              py=py+dr
              pz=pz+dr
           case(3)
              px=px-dr
              py=py-dr
              pz=pz+dr
           case(4)
              px=px+dr
              py=py-dr
              pz=pz+dr
           case(5)
              px=px+dr
              py=py+dr
              pz=pz-dr
           case(6)
              px=px-dr
              py=py+dr
              pz=pz-dr
           case(7)
              px=px-dr
              py=py-dr
              pz=pz-dr
           case(8)
              px=px+dr
              py=py-dr
              pz=pz-dr
           end select
        case(4)
           ! determine new position
           call random_number(rn)
           rn=rn*2.0*pi
           call random_number(rn2)
           rn2=rn2*2.0*pi
           d1=cos(rn2)
           px=px+dist*cos(rn)*d1
           py=py+dist*sin(rn)*d1
           pz=pz+dist*sin(rn2)
        end select
     end if
     return
  case(2)
     ! straight line

     ! determine the random position between
     ! (x1,y1,z1) and (x2,y2,z2)
     call random_number(rn)
     px=geom_obj(gid)%x1+&
          (geom_obj(gid)%x2-geom_obj(gid)%x1)*rn
     py=geom_obj(gid)%y1+&
          (geom_obj(gid)%y2-geom_obj(gid)%y1)*rn
     pz=geom_obj(gid)%z1+&
          (geom_obj(gid)%z2-geom_obj(gid)%z1)*rn
     
     ! determine new position in isotropic emission
     ! This must be fixed later.
     call random_number(rn)
     rn=rn*2.0*pi
     call random_number(rn2)
     rn2=rn2*2.0*pi
     d1=cos(rn2)
     px=px+dist*cos(rn)*d1
     py=py+dist*sin(rn)*d1
     pz=pz+dist*sin(rn2)
  case(3)
     ! x=a,left
     px=geom_obj(gid)%x1+dist
     call random_number(rn)
     py=rn*(geom_obj(gid)%y2-geom_obj(gid)%y1)+geom_obj(gid)%y1
     call random_number(rn)
     pz=rn*(geom_obj(gid)%z2-geom_obj(gid)%z1)+geom_obj(gid)%z1
  case(4)
     ! x=a,right
     px=geom_obj(gid)%x1-dist
     call random_number(rn)
     py=rn*(geom_obj(gid)%y2-geom_obj(gid)%y1)+geom_obj(gid)%y1
     call random_number(rn)
     pz=rn*(geom_obj(gid)%z2-geom_obj(gid)%z1)+geom_obj(gid)%z1
  case(5)
     ! y=b,left
     py=geom_obj(gid)%y1+dist
     call random_number(rn)
     px=rn*(geom_obj(gid)%x2-geom_obj(gid)%x1)+geom_obj(gid)%x1
     call random_number(rn)
     pz=rn*(geom_obj(gid)%z2-geom_obj(gid)%z1)+geom_obj(gid)%z1
  case(6)
     ! y=b,right
     py=geom_obj(gid)%y1-dist
     call random_number(rn)
     px=rn*(geom_obj(gid)%x2-geom_obj(gid)%x1)+geom_obj(gid)%x1
     call random_number(rn)
     pz=rn*(geom_obj(gid)%z2-geom_obj(gid)%z1)+geom_obj(gid)%z1
  case(7)
     ! z=c,left
     pz=geom_obj(gid)%z1+dist
     call random_number(rn)
     px=rn*(geom_obj(gid)%x2-geom_obj(gid)%x1)+geom_obj(gid)%x1
     call random_number(rn)
     py=rn*(geom_obj(gid)%y2-geom_obj(gid)%y1)+geom_obj(gid)%y1
  case(8)
     ! z=c,right
     pz=geom_obj(gid)%z1-dist
     call random_number(rn)
     px=rn*(geom_obj(gid)%x2-geom_obj(gid)%x1)+geom_obj(gid)%x1
     call random_number(rn)
     py=rn*(geom_obj(gid)%y2-geom_obj(gid)%y1)+geom_obj(gid)%y1
  end select

  return

end subroutine set_source_pos
